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Differential  Geometric  Tools  for  Image  Sensor  Fusion 

Final  Report 

AFOSR  STTR  Phase  I  contract  #  F49620-98-C-0047 


Diego  A,  Socolinsky,  Lawrence  B.  Wolff 


Equinox  Corporation 
9  West  57th  Street 
New  York,  New  York  10019 

Abstract 

This  report  describes  work  performed  under  AFOSR  STTR  Phase  I  contract  #  F49620-9S- 
C-0047  from  September  1,  1998  through  August  31,  1999.  We  present  a  new  formalism  for  the 
treatment  and  understanding  of  multispectral  images  and  multisensor  imagery  based  on  first 
order  contrast  information.  Although  little  attention  has  been  paid  to  the  utility  of  multispectral 
contrast,  we  develop  a  theory  for  multispectral  contrast  that  enables  us  to  produce  an  optimal 
grayscale  \isualization  of  the  first  order  contrast  of  an  image  with  an  arbitrary  number  of  bands. 
We  demonstrate  how  our  technique  can  reveal  significantly  more  interpretive  information  to  an 
image  analyst,  who  can  use  it  in  a  number  of  image  understanding  algorithms.  Existing  grayscale 
visualization  strategies  are  reviewed.  A  variety  of  experimental  results  are  presented  to  support 
the  performance  of  the  new  method. 


L  Introduction 

The  advent  of  new  remote  sensing  and  imaging  technologies  provides  us  with  ever  increasing  vol- 
imes  of  multispectral  data.  Faced  with  this  information  explosion,  it  has  become  necessary  to 
levelop  methods  for  analysis  of  such  high  dimensional  datasets.  One  key  aspect  of  this  process  is 
;he  visualization  of  multispectral  data,  to  be  used  for  photointerpretation.  This  allows  an  image 
malyst  to  determine  regions  of  interest  and  important  features  in  the  image  for  further  analysis 
)r  segmentation.  In  order  to  take  full  advantage  of  the  human  visual  system,  a  Red-Green-Blue 
composite  image  is  usually  generated  from  the  data  by  one  of  a  number  of  statistical  methods  which 
VC  review  in  section  2.  We  aim  to  produce  a  one-band,  grayscale  visualization  image  from  a  given 
nultispectral  dataset.  We  propose  to  do  this  in  such  a  way  as  to  preserve  as  much  local  image 
:ontrast  ‘feature  information’  as  possible. 

Computation  of  contrast,  which  includes  computation  of  gradient  and  zero-crossings,  has  been 
ised  in  computer  vision  as  one  of  the  primary  methods  for  extracting  grayscale  and  color  features 
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16,  2,  4,  10,  6].  It  seems  plausible,  therefore,  that  the  correct  way  to  compare  versions  of  the 
same  image  in  terms  of  feature  information  is  through  their  contrast  content.  Contrast,  however,  is 
leSned  a  priori  only  for  grayscale  images,  so  we  cannot  readily  compare  multiband  images  amongst 
shemselves  or  multiband  images  with  grayscale  images.  The  first  step  is  thus  to  define  contrast  for 
i  multiband  image.  This  is  achieved  in  section  3.1  through  the  introduction  of  a  differential  form 
m  the  image,  computed  in  terms  of  the  spectral  map  and  a  metric  defined  on  photometric  space. 
This  reduces  to  the  standard  notion  of  contrast  in  grayscale  images. 


Once  contrast  has  been  defined  for  an  arbitrary  image,  it  is  natural  to  ask  which  grayscale  image 
nost  closely  matches  the  contrast  information  of  a  given  multiband  image.  Or,  how  should  we 
:onvert  a  multiband  image  to  grayscale  while  preserving  as  much  contrast  information  as  possible? 
This  problem  is  taken  up  in  section  4,  where  we  show  the  mathematical  formulation  of  the  problem, 
logether  with  its  solution  and  experimental  results. 


It  should  be  noted  that  the  solution  to  this  problem  has  multiple  applications.  For  the  remote 
sensing  community,  this  algorithm  provides  a  visualization  tool  for  realizing  the  full  edge  infor- 
nation  content  in  raultispectral  images,  such  as  those  obtained  through  satellite  imaging.  Such 
ligh-dimensional  photometric  data  is  not  easily  tractable  by  traditional  methods;  in  this  context 
)ur  reduction  method  yields  a  useful  data  analysis  tool.  In  medical  imaging,  image  fusion  can  be 
ised  to  simultaneously  visualize  multiple  data  modalities  such  as  CT,  PET  and  MRI  [19]. 


2  Review  of  existing  techniques 

^erhaps  the  simplest  possible  transformation  from  a  multispectral  image  to  a  grayscale  image  is 
Lveraging  of  the  spectral  bands.  This  produces  a  visualizable  image  which  contains  information 
rom  all  the  bands  in  a  unified  way.  However,  as  is  easily  seen,  this  method  fails  to  take  into  account 
<  my  measure  of  the  information  content  in  the  dataset.  A  minor  modification  can  be  obtained  by 
i  onsidering  a  weighted  aviarage,  where  different  bands  will  contribute  differently  to  the  final  result, 

( lepending  on  some  pre-aseigned  assesment  of  their  relative  relevance  in  the  overall  image.  Since  it 
:  nay  be  dilSicult  or  even  impossible  to  determine  a  priori  which  bands  should  be  emphasized  over 
( others,  this  method  suffers  from  similar  problems  as  unweighted  averaging.  Furthermore,  relative 
j  ;rayvalues  in  different  bands  may  be  such  that  features  are  completely  obliterated  by  this  process, 
:  or  example  consider  averaging  two  black  and  white  checkerboard  patterns  with  grayvalues  reversed. 

In  order  to  overcome  the  shortcomings  of  averaging  methods,  we  can  take  into  account  statistical 

3  nformation  about  the  multispectral  image.  Principal  Component  Analysis  (PCA)  achieves  this  by 
i  onsidering  an  n-band  image  as  a  set  of  vectors  in  an  n-dimensional  vector  space,  A  grayscale 
^  isualization  is  obtained  from  a  multipectral  image  by  projecting  the  entire  distribution  of  spectral 
'  alues  onto  the  line  spanned  by  the  eigenspace  of  the  covariance  matrix  with  largest  eigenvalue, 
3  nd  then  perhaps  re-scaling  the  result  to  fit  the  dynamic  range  of  the  output  device  (printer, 
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Figure  1:  Angle  between  diagonals  and  coordinate  axes  as  a  function  of  spectral  dimension. 


nonitor,  etc).  This  technique  has  been  extensively  used  for  visualization  of  remote  sensing  imagery 
13],  as  well  as  multimodal  medical  imaging  [18].  The  difference,  between  PC  A  and  weighted 
iveraging  is  that  the  line  onto  which  we  project  is  chosen  ahead  of  time  in  the  latter,  whereas 
n  the  former  it  is  determined  by  the  global  statistics  of  the  particular  image  at  hand.  However 
lince  both  methods  have  a  common  geometric  foundation,  they  share  a  common  problem.  To  see 
;his  clearly  we  resort  to  the  following  argument  [12].  It  is  easy  to  see  (Figure  1)  that  the  cosine 
)f  the  angle  6  between  any  diagonal  vector  in  an  n-dimensional  vector  space  and  any  one  of  the 
:oordinate  axis  is  given  by  co8(0)  =  Hence  as  the  dimension  increases,  diagonals  tend  to  become 
orthogonal  to  the  coordinate  axes.  It  follows  that  upon  projecting  the  spectral  measurements 
onto  a  fixed  axis  or  a  principal  axis  in  photometric  space,  the  contrast  between  adjacent  pixels  is 
i  always  foreshortened,  and  it  follows  from  the  previous  remark  that  this  phenomenon  becomes  more 
1  evere  as  the  dimensionality  of  the  data  increases.  Any  method  based  on  linear  projections  will  be 
I  Aversely  affected  by  this  situation. 

Multiresolution  methods  based  on  pyramidal  decompositions  [22,  3,  21]  and  wavelet  transforms 
9]  have  been  reported.  A  common  feature  to  all  these  methods  is  a  selection  rule  which  determines 
^  /hich  band  of  the  multiband  image  is  'most  relevant’  in  a  neighborhood  of  a  given  pixel;  the  features 
i  f  the  selected  band  are  then  incorporated  into  the  fused  image  through  various  procedures.  By 
1  Lsing  such  a  selection  rule,  these  methods  implicitly. assume  that  there  is  only  one  dominant  band  at 
i  ach  pixel.  It  is  well-known  that  multispectral  imagery  often  exhibits  large  inter-band  correlation, 
f  o  this  assumption  is  often  violated.  The  consequence  is  that  such  methods  do  not  allow  for  small 
1  tut  consistent  contrast  features  across  different  bands  to  compound  and  form  more  salient  features 
i  Cl  the  fused  image.  The  loss  of  potential  contrast  resulting  from  such  a  selection  rule  is  easily  seen 
t  D  be  proportional  to  the  square  root  of  the  number  of  bands. 

Recently  a  number  of  other  visualization  methods  have  been  proposed,  most  notably  those 
based  on  self-organizing  maps  [l]  and  optimal  projection  maps  [11].  While  very  different  in  their 
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ideas  and  implementation,  both  methods  rely  on  global  image  statistics  to  produce  a  linear  or 
aon-linear  projection  of  multispectral  apace  onto  a  grayscale  axis.  The  use  of  global  statistics  in 
visualization  technique  has  undesireable  side-effects,  which  we  exemplify  in  Section  4.3. 


3  Contrast  in  a  multiband  image 

\  standard  reference  for  all  geometric  notions  used  in  this  section  is  [20]. 

3.1  The  contrast  form 

3y  a  multiband  image,  we  mean  a  rectangle  C  together  with  a  spectral  map  s  :  fl  — >  P",  where 
P”  denotes  rt-dimensional  photometric  space.  We  will  assume  that  P"  is  given  an  ‘extended’  RGB 
coordinate  system,  in  which  each  band  takes  on  values  between  0  and  M  <  cc.  The  extension  to 
)ther  coordinate  systems  is  straightforward.  In  this  context,  a  grayscale  image  is  a  map  a  ;  fl  — >  P^. 
We  let  P”  have  an  arbitrary  Riemannian  metric  g,  which  can  be  used  to  address  the  issue  of  sensor 
loise  as  in  [23,  7],  to  introduce  psychophysically  correct  metrics  [8],  or  simply  to  manually  influence 
;he  fusion  (see  example  in  Section  4.3).  We  expand  on  this  in  section  4.2.  Finally,  if  the  different 
>ands  in  the  image  are  the  result  of  measurements  taken  with  different  sensors  and/or  at  different 
limes,  then  we  assume  that  these  images  have  been  properly  registered  with  respect  to  eadi  other 
ind  brought  to  a  common  spatial  resolution. 

For  the  remainder,  let  s  :  fl  — >  P’^  be  a  multiband  image,  let  p  be  a  point  in  fl,  and  v  an 
ixbitrary  unit  vector  in  .  In  analogy  with  the  grayscale  case,  we  seek  to  define  the  contrast  in 
!  at  p  in  the  direction  of  v  as  the  speed  of  spectral  variation  within  the  image  in  that  direction. 
Do  do  so,  consider  the  following  construction.  Let  7  ;  [~c,  e]  -+  0  be  a  curve  defined  on  a  small 
nterval,  such  that  7(0)  =  p  and  7^(0)  =  v.  The  speed  of  spectral  variation  at  p  in  the  direction 
>f  V  is  given  by  the  magnitude  of  the  vector  s*(i')  =  o  7) (t)  1^=0.  as  evaluated  by  the  metric 
m  P”.  Note  that  since  the  metric  of  P"  is  not  assumed  to  be  constant,  this  magnitude  must  be 
jvaluated  with  respect  to  the  metric  at  3(p}.  In  the  language  of  differential  geometry,  the  vector 
'*(v)  is  called  the  pushforward  of  v  by  s,  and  its  expression  in  local  coordinates  is  given  by 

s*(u)  =  (1) 

vhere  Jp  is  the  Jacobian  matrix  of  s  at  the  point  p.  Let  g,  denote  the  matrix  for  5  at  g  G  P”  with 
aspect  to  a  coordinate  system.  Then  the  contrast  of  s  at  p  in  the  direction  of  v  is  given  by  the 
luantity 


{Jpv)  gs(p){Jpv)  —  V  {Jpgs(p)Jp)v. 


(2) 
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Let  5*  -  ( Jp).  It  follows  from  Equation  (1)  that  in  a  coordinate  system  the  components  of 
ffp  axe  given  by 


(ffp)u  -  S  (3) 

Finally,  define  x^ip)  =  9p  to  be  the  image  contrast  form  of  s.  Thus  the  contrast  of  s  at  p  e  n  in 
the  direction  of  w  is  given  by  x^(p)(u).  Once  again,  in  the  language  of  differential  geometry,  the 
differential  form  x  “  called  the  pullback  of  g  by  s.  In  a  coordinate  system,  the  image  contrast  form 
can  be  expressed  as 


(4) 


which  upon  evaluation  on  a  vector  v  =  yields  the  non-negative  real  number 


x=(P)w  =  Etfe(P.)«i|^“V. 

i,i=i  t,<=i 


(5) 


This  differential  form  encodes  all  first  order  contrast  information  about  a  multispectral  image.  Any 
false  color  image  visualization  whose  corresponding  image  contrast  differential  form  is  identical  to 
that  for  a  multispectral  image,  replicates  its  first  order  contrast  information. 

Note  that  if  all  bands  are  brought  to  a  common  resolution  prior  to  the  computation  of  the 
contrast  form  (by  re-sampling,  interpolation  or  any  other  suitable  method),  then  the  choice  of  g  is 
independent  of  the  original  resolution  since  p  is  a  metric  on  photometric  space,  not  on  the  image 
plane.  However,  if  one  chooses  to  carry  the  computations  in  the  original  resolutions,  then  g  can 
be  used  to  compensate  for  the  disparity,  as  long  as  the  difference  in  resolution  is  uniform  in  both 
jpatial  coordinates.  Thus,  for  example  if  one  band  has  half  the  resolution  another,  then  we  can 
ise  a  metric  that  assigns  the  lower  resolution  band  half  as  much  weight  as  to  the  other.  It  is  more 
Jtraightforward  and  efficient  to  bring  all  bands  to  a  common  resolution  prior  to  computing  the 
contrast  form,  and  so  we  use  this  method  in  general- 

Prom  Equation  (3)  we  see  that  x^(p)  is  a  symmetric  matrix  with  real  entries,  therefore  its 
agenvalues  axe  both  real  and  non-negative.  Let  Xp  denote  the  largest  eigenvalue  of  x^ip)-  We 
lefine  the  absolute  contrast  of  s  at  p  G  fl  to  be  equal  to  and  we  say  that  the  eigenspace 
X  (p)  corresponding  to  Xp  is  the  direction  of  maximal  contrast  at  p.  Note  that  this  direction 
loes  not  have  a  preferred  orientation  a  priori,  i.e.  the  eigenspace  is  a  line  without  a  preassigned 
>rientation. 


(Let  us  pause  for  a  moment  and  see  how  these  definitions  simplify  in  the  case  where  fl  is  a 
rayscale  image  and  is  given  the  standard  Euclidean  metric.  Using  coordinates  x.  y  on  fl  and 
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l6tting  subscripts  denote  partial  differentiation,  we  see  that  the  Jacobian  matrix  of  s  is  Sj,) , 
and  therefore 


We  can  readily  compute  the  eigenvalues  of  this  matrix  to  be  0  and  |Vsp.  Thus  Ap  =  |V5(p)|,  and 
the  direction  of  maximal  contrast  is  the  line  spanned  by  Vs(p).  Hence  we  recover  the  standard 
notion  of  grayscale  contrast,  modulo  orientation.  In  this  case,  note  that  the  direction  of  maximal 
contrast  can  be  given  the  orientation  determined  by  Va.  For  a  general  multiband  image,  there 
is  no  canonical  orientation  for  the  direction  of  maximal  contrast;  this  constitutes  one  of  the  main 
differences  between  multiband  and  grayscale  images  in  terms  of  contrast. 

It  follows  from  the  previous  considerations  that  two  images  cairry  the  same  contrast  information 
if  and  only  if  they  have  the  same  contrast  form.  Note  that  the  contrast  form  of  a  multiband  image 
will,  in  general,  have  two  non-zero  eigenvalues.  Since  the  contrast  form  of  a  grayscale  Image,  shown 
above,  always  has  one  null  eigenvalue,  we  have  that  the  contrast  information  in  a  multiband  image 

cannot,  in  general,  be  exactly  reproduced  by  a  grayscale  image.  This  observation  will  be  necessary 
in  section  4.1. 

3.2  The  contrast  vector  field 


The  definition  of  contrast  given  in  the  previous  section  is  not  quite  sufficient  to  tackle  the  optimal 
(visualization  problem.  This  is  so  precisely  because,  as  we  noted  before,  there  is  no  canonical 
orientation  for  the  direction  of  maximal  contrast  at  a  given  point  in  a  multiband  image.  We  must 
remedy  this  by  introducing  such  an  orientation  in  a  well-defined  manner. 

Let  X(p)  *=  dtst(0,  s(p)),  where  0  is  the  point  in  P”  all  of  whose  coordinates  are  zero,  and  dtsi  is 
;he  geodesic  distance  function  for  the  metric  of  P".  The  function  2  :  -o  M  is  the  spectral  intensity 
unction.  Note  that  in  the  case  of  Euclidean  photometric  space  we  simply  have  2(p)  =  3«(p)^- 

:n  general,  2  induces  an  ordering  on  P",  given  by  x  >  y  if  2(ir)  >  2(y),  for  x,y  e  P.  The  function 
Hat  could  in  principle  be  replaced  here  by  any  non-negative  function  :  P”  -4  R  whose  level  sets 
tre  hypersurfacee  foliating  F*,  satisfying  the  relation  ^6-^([0,o])  C  ^“^[0,6]),  for  all  0  <  c  <  b. 
The  spectral  intensity  defined  using  ^  instead  of  dist  would  thus  induce  a  different  ordering  of  the 
let  of  colors. 

Now,  we  construct  the  contrast  vector  field  V  as  follows.  Let  Vp  be  the  unique  vector  of  length 
along  the  direction  of  maximal  contrast  at  p,  for  which  Vvj,2  >0.  If  P"  is  given  an  arbitrary 
liemannian  metric,  this  procedure  can  be  achieved  by  considering  the  sign  of  the  Riemannian 
:nner  product  between  Vp  and  the  outward  unit  normal  to  the  geodesic  sphere  of  radius  2(p)  at 
.  (p)  instead  of  the  sign  of  (Vtt,2)(p).  If  a  different  function  <f>  is  chosen,  then  ‘geodesic  sphere’ 
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figure  2:  The  direction  of  the  contrast  vector  at  p,  is  the  image  direction  at  the  point  p  along  the 
argest  eigenvector  of  that  produces  a  positive  variation  in  I. 

I  hould  be  replaced  by  ‘level  set  of  in  the  previous  sentence.  Essentially,  we  orient  V  to  point  in 
‘  he  direction  of  increasing  spectral  intensity. 

It  is  quite  important  to  allow  for  different  choices  of  4>,  as  can  clearly  be  seen  by  considering  an 
( ixample.  If  we  take  a  2-band  image  whose  bands  are  very  negatively  correlated,  the  orientation  of 
will  be  inconsistent  along  image  edges.  An  artificial  situation  exhibiting  this  behavior  can  be 
i  onstructed  by  taking  the  two  bands  to  be  Z+t?  and  tj',  where  /  is  any  intensity  function,  and  t) 
\  nd  rf  are  random  gaussian  fields  representing  image  noise.  In  this  case,  the  orientation  of  Vp  using 
1  he  ordering  induced  by  the  Euclidean  metric  is  entirely  dependent  on  the  noise,  and  completely 
unrelated  to  the  image  features.  One  valid  choice  of  <!>  in  this  case  is  any  of  the  projections 
( onstructed  in  [11]  (but  note  that  these  projections  may  not  preserve  the  ‘natural’  relation  of  dark 
j  nd  bright  in  the  image). 

The  contrast  vector  field  V  constructed  above,  encodes  the  first  order  spectral  contrast  infer- 
]  aation  of  Q,  together  with  intensity  information.  It  constitutes  a  bridge  between  the  multispectral 
<  nd  grayscale  realms. 

- 1  Grayscale  visualization  of  local  contrast 

]  'or  the  analytic  aspects  of  this  section,  the  reader  may  refer  to  [5].  The  numerical  methods  can  be 
ljund  in  [14]. 
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4.1  Mathematical  formulation 

In  light  of  the  definitions  of  section  3.2,  the  statement  of  the  problem  is  now  the  following:  Given  a 
multiband  image  a  :  Q  — P”  with  contrast  vector  field  V,  find  the  grayscale  image  whose  contrast 

vector  field  is  closest  to  V.  Equivalently,  we  seek  the  function  /  :  fl  K  whose  gradient  is  closest 
to  V. 

To  find  such  a  function,  we  would  like  to  solve  the  equation  V/  =  V.  However,  this  equation 
will  in  general  have  no  solution,  since  it  follows  from  the  remarks  at  the  end  of  section  3.1  that  V 
need  not  be  integrable.  This  means  that  in  general  it  will  not  be  possible  for  a  grayscale  image  to 
exactly  reproduce  the  contrast  information  of  a  multispectral  image.  We  propose  instead  to  find 
the  function  /  for  which  the  following  functional  is  minimized 

I  J\'^f-V\^dzdy.  (7) 

n 

The  Euler-Lagrange  equation  for  this  functional  can  be  easily  found  to  be 

f  A/  =  divV,  on  U. 

^  V/  •  n  =  y  -  n,  on  30, 

where  is  the  outward  unit  normal  to  O.  There  is  no  natural  way  to  assign  Dirichlet  boundary 
conditions,  so  we  choose  Neumann  conditions,  consistent  with  the  idea  that  V  should  approximate 
bhe  gradient  of  the  solution  image  intensity  function. 

We  would  like  to  point  out  that  this  formulation  represents  a  new  paradigm  in  multispectral 
mage  visualization.  The  standard  approach  to  the  problem  has  been  to  seek  a  projection  (usually 
inear)  from  photometric  space  onto  a  one-dimensional  grayscale  axis.  Our  formulation  omits  the 
projection  and  instead  seeks  the  best  grayscale  image  itself.  The  main  consequence  of  this  is  that 
ve  have:  much  more  freedom  to  reproduce  contrast  variations.  On  the  other  hand,  a  side  effect  is 
;hat  pixels  with  the  same  photometric  values  in  the  multispectral  image  may  not  have  the  same 
grayscale  value  in  our  visualization.  We  have  not  found  this  to  be  a  problem,  since  it  simply  reflects 
ihe  fact  that  in  our  formulation,  contrast  is  a  local  quantity.  In  the  conclusion  we  outline  a  solution 
hat  produces  a  unique  mapping  of  grayvalues  at  the  expense  of  reduced  contrast  fidelity.  Also, 
n  section  4.3  we  show  how  one  can  combine  the  utility  of  the  new  method  with  some  desireable 
properties  of  statistical  projections,  to  obtain  high-contrast  color  composites  whose  color  scheme 
irovides  a  unique  correspondence  between  raw  photometric  values  and  image  colors. 

There  are  a  number  of  methods  for  the  solution  of  Poisson’s  Equation  (8)  with  Neumann 
loundaiy  conditions  on  a  rectangle.  A  simple  iterative  scheme  based  on  the  standard  five-point 
1  ipproximation  for  the  Laplacian  can  be  written  for  a  discrete  pixel  grid  of  dimensions  [0,  i]  x  [0,  J]. 
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Letting  /°  be  any  initial  guess,  we  write 

A'j‘ = 4 + i  .  (9) 

for  (i,J)  €  [0,7]  X  [0,  J].  For  interior  pixels  (i,j)  €  [1,/  —  1]  x  [1,  J  —  1],  we  use 

=  fUij  +  fi-U  +  fij^i  +  -  ^flj-  ( 10) 

The  divergence  term  must  be  discretized  appropriately  to  avoid  artifacts  in  the  solution  image. 
For  the  Laplacian  operator  in  (10),  we  construct  the  contrast  vector  field  using  forward  difference 
derivatives  in  the  computation  of  The  divergence  is  then  computed  using  backward  difference 
derivatives  of  the  components  of  V.  In  this  way  we  insure  that  the  correct  discrete  solution  is  found 
in  the  case  where  V  is  itself  the  gradient  of  some  function. 

At  the  boundary  pixels,  Equation  (10)  is  undefined.  Here,  the  Neumann  condition  in  (8)  is  used 
to  modify  (10).  For  example,  if  (t,  J),  1  <  t  <  7  -  1,  is  a  pixel  on  the  lower  horizontal  boundary 
segment  of  the  image,  the  boundary  condition  becomes  fi,j+i  -  fij  =  Vi^j,  which  impUes  the 
Lapla/cian  operator 


^fi.J  —  +  fi-l,J  +  ftj-1  +  ~ 


(11) 


The  boundary  condition  can  similarly  be  used  to  modify  the  computation  of  V  and  div  V  at  bound¬ 
ary  pixels.  Note  that  convergence  is  guaranteed  regardless  of  the  initial  guess,  and  all  solutions 
agree  upto  an  overall  additive  constant. 

A  common  modification  of  (9)  in  the  interest  of  computational  speed  is  the  addition  of  over-  or 
under-relaxation  to  speed  convergence  [14,  15].  In  this  case  (9)  becomes 


=(i— 


(12) 


for  some  choice  of  relaxation  parameter  0  <  w  <  2.  It  can  be  shown  [14]  that  for  the  discretization 
above  on  the  pixel  grid  [0, 7]  x  [0,  J],  the  ideal  relaxation  parameter  is  given  by 


,  cos(7r/7)  +  cos(7r/J) 

for  p  = - i ^ 


Since  we  normally  work  on  large  pixels  grids  with  7,  J  >  100,  the  relaxation  parameter  is  very  close 
to  2,  and  large  speed  improvement  results  from  over-relaxation. 
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4.2  Noisy  data 

[t  may  be  the  case  that  all  bands  of  a  multispectral  image  have  comparable  resolutions,  but  still 
each  band  has  different  noise  characteristics.  This  may  be  due  to  a  number  of  reasons;  different 
technologies  are  used  for  sampling  at  different  wavelengths,  penetration  of  the  atmosphere  at  dif¬ 
ferent  wavelengths  is  not  uniform,  cloud  cover  in  a  satellite  image  may  introduce  noise  only  in  the 
visible  range,  but  not  in  the  infraxed,  etc.  In  some  of  these  cases,  noise  characteristics  of  the  sensor 
at  each  band  may  be  known  a  priori.  The  approach  described  in  [23]  allows  us  to  temper  the  effect 
of  noise  in  this  situation:  the  n  x  n  covariance  matrix  associated  to  the  sensor  is  computed  (or 
interpolated)  at  each  point  in  n-dimensional  photometric  space,  yielding  a  positive  deffnite  bilinear 
form  which  can  be  thought  of  as  a  Riemannian  metric  on  photometric  space.  This  metric  encodes 
the  noise  characteristics  of  the  sensor  in  such  a  way  that  unit  (in  the  Euclidean  norm)  tangent 
vectors  in  the  direction  of  noisier  bands  will  have  shorter  lengths  than  those  in  the  direction  of 
bands  of  higher  fidelity.  We  use  this  as  the  metric  g  in  Equation  (2)  and  proceed  with  the  rest 
of  the  algorithm  with  no  further  changes.  The  overall  effect  is  to  have  the  higher  fidelity  bands 
contribute  more  to  the  grayscale  composite  than  those  which  are  known  to  be  noisier. 

If  we  have  no  o  priori  knowledge  of  the  sensor  noise  characteristics,  but  we  assume  independent 
identically  distributed  noise  in  all  bands,  with  zero  mean  and  distribution  symmetric  about  the 
mean,  then  we  can  use  a  hybrid  of  the  method  above,  and  band  averaging.  Let  be  a  2- 
dimensional  Gaussian  distribution  with  zero  mean  and  standard  deviation  tr .  Compute  the  low- 
and  high-frequency  components  of  the  multiband  image  s  as 

Itf  =  Ga  *  s,  arid  h(,  =  s  —  Gg  *  s.  (13) 

Let  Xff  til®  contrast  form  of  the  low-pass  filtered  image  1^^,  and  fa  be  the  grayscale  fused 
visualization  resulting  from  th®  method  introduced  above.  Lastly,  construct  the  final 

visualization  as 

+  (14) 

The  assumptions  on  the  noise  distribution  are  such  that  high  frequency  noise  in  the  first  term  of 
(14)  will  tend  to  cancel  out,  while  lower  frequency  features  from  each  band  will  combine  within  the 
second  term  to  produce  a  fusion  of  the  low-pass  filtered  bands.  The  parameter  a  clearly  controls 
the  frequency  threshold  past  which  we  consider  the  data  to  be  too  noisy,  and  in  the  limit  as  a  0 
we  recover  the  original  procedure,  with  no  noise  attenuation. 

A  more  involved  procedure  based  on  an  adaptive  metric  is  developed  in  [7],  A  complete  descrip¬ 
tion  of  the  method  is  outside  the  scope  of  the  current  article,  but  we  include  an  example  image  here 
to  deinonstrate  the  effect  of  the  metric  in  the  context  of  noise  suppression.  The  images  in  Figure  3 
were  obtained  through  the  application  of  the  technique  introduced  in  this  paper,  to  a  4-band  aerial 
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image  with  rather  severe  noise  corruption.  The  image  in  3(a)  was  created  using  a  Euclidean  metric 
whereas  the  one  in  3(b)  exploits  the  adaptive  construction  in  [7]. 


(a)  (b) 


Figure  3:  Fusion  of  4-band  aerial  image  corrupted  by  noise:  (a)  with  Euclidean  metric,  (b)  with 
adaptive  metric  from  [7] 


4.3  Experimental  results 

We  performed  a  series  of  experiments  to  validate  our  theoretical  approach,  using  7-dimensional 
data  from  the  EOS  AT  thematic  mapper  and  12-dimensional  remote  sensing  imagery,  and  2 10- band 
data  from  the  HYDICE  project,  as  well  as  a  variety  of  multimodal  medical  images.  All  the  fused 
images  shown  in  this  section  were  computed  using  a  Euclidean  metric,  unless  otherwise  noted. 

Figures  4(a)  and  4(b)  show  the  result  of  applying  PCA  and  our  algorithm,  respectively,  to  a 
12-band  image.  A  number  of  corresponding  regions  have  been  selected  in  both  images  to  highlight 
some  of  the  differences  between  them.  Note  how  new  features  axe  visible  in  our  visualization  in 
the  areas  labeled  1,  2.  The  continuation  of  a  road  which  gets  lost  in  the  PCA  image  is  visible 
in  area  4.  More  contrast  detail  is  present  in  area  3,  to  the  extent  that  it  is  possible  to  identify 
a  distinctly  rectangular  region  almost  invisible  in  the  PCA  visualization.  Area  5  shows  higher 
resolution  of  objects  in  the  image.  The  clear  advantage  of  our  visualization  technique  would  allow 
an  image  analyst  to  more  accurately  and  reliably  select  regions  of  interest  in  the  image  for  further 
processing  by  image  understanding  algorithms.  Since  these  algorithms  become  significantly  more 
time  consuming  as  the  number  of  bands  grows,  the  ability  to  quickly  reduce  the  axeas  to  which 
they  must  be  applied  results  in  increased  processing  efficiency. 
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Figures  5  and  6  show  more  detailed  views  of  our  methodology  at  work  versus  principal  com¬ 
ponent  analysis.  We  can  see  that  not  only  does  our  technique  highlight  more  detail,  but  it  even 
incovers  features  that  were  obliterated  by  the  statistical  method.  It  is  hard  to  overestimate  the 
idvantage  this  affords  the  image  analyst.  In  Figure  5  we  see  how  a  number  of  small  islands  along 
the  coastline  of  Lake  Worth,  Florida,  dissapear  from  the  statistical  visualization,  while  they  are 
easily  visible  in  ours.  Likewise  Figure  6  shows  how  we  recover  more  of  the  original  man-made 
(probably  airport  runways)  and  natural  structures  from  a  9  band  EOSAT  image  than  PCA. 

Figures  7(a)  and  7(b)  show  the  results  of  applying  PCA  and  oux  algorithm  to  an  artificial  movie 
sequence  produced  by  sliding  a  window  over  area  2  of  the  12-band  image  in  Figure  4,  respectively. 
First  note  the  greatly  enhanced  level  of  feature  information  conveyed  by  the  images  obtained  by 
Dur  method.  In  fact  the  roads  and  divisions  between  land  plots  axe  hardly  visible  in  7(a).  We 
also  see  that  the  same  physical  feature  looks  different  in  different  frames  of  the  sequence  in  Figure 
7(a).  The  grayscale  values  for  a  given  feature  do  not  remain  constant  through  the  sequence,  as 
they  depend  on  the  global  statistics  of  the  particular  image  for  which  the  principal  components 
were  computed.  We  should  mention  that  this  shortcoming  is  not  unique  to  PCA;  it  affects  any 
visualization  method  which  exploits  the  global  statistics  of  the  image.  In  sharp  contrast  to  this 
phenomenon,  the  movie  sequence  produced  by  our  visualization  algorithm  is  consistent  in  terms  of 
grayscale  values.  This  property  makes  it  possible  for  us  to  create  consistent  optimal  visualizations 
of  multispectral  video.  We  are  not  aware  of  any  other  methods  which  axe  currently  capable  of  this. 

Grayscale  consistency  on  overlapping  regions  is  an  asset  for  multispectral  image  registration 
as  well.  Many  manual  and  semi-automatic  registration  methods  rely  on  an  operator  selecting  a 
number  of  matching  features  from  both  images.  These  features  are  then  used  as  tie-points  for  a 
discrete  optimization  algorithm.  When  distinctive  features  occur  across  different  spectral  bands,  it 
is  advantageous  for  the  operator  to  use  band  composites  for  feature  selection.  However,  if  statistical 
methods  are  used  to  create  the  composites,  then  the  overlapping  region  between  the  images  may 
appear  quite  different  in  each  image.  This  renders  the  tie-point  selection  process  more  difficult 
and  less  reliable.  By  using  our  visualization  algorithm  it  is  possible  to  obtain  images  with  both 
consistent  shading  in  the  overlapping  region  and  rich  feature  content. 

Figure  8  shows  the  relative  performance  of  our  method  and  the  multiresolution  wavelet  method 
in  [9]  for  the  fusion  of  chest  CT.  We  used  two  contrast  windows  of  body  CT  to  produce  these 
unified  visualizations.  Note  how  the  wavelet  method  yields  a  lower  contrast  image,  in  which  many 
of  the  details  appear  washed  out,  while  our  method  produces  crisper  detail  throughout.  No  post¬ 
processing  was  done  on  either  image. 

Let  us  see  by  example  how  particular  choices  of  metric  in  photometric  space  can  be  used  to 
achieve  various  desired  effects.  Figure  9  shows  three  modalities  of  brain  MRI  of  the  same  patient. 
In  Figure  10(a)  we  have  the  result  of  using  a  Euclidean  metric  in  Equation  (5)  for  the  fusion  of 
these  three  modalities-  Note  among  other  details,  how  the  tumor  from  9(b)  and  the  skull  from 
9(a,c)  are  both  clearly  visible.  To  create  the  image  in  10(b),  we  used  a  metric  which  emphasizes 
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contrast  in  the  T2  weighted  MRI  for  dark  pixels  (10-50  gray  counts  out  of  256),  which  has  the 
effect  df  increasing  contrast  between  gray  and  white  matter  in  the  fused  image  without  affecting,  for 
example,  the  visibility  of  the  tumorous  lesion.  The  opposite  effect  is  displayed  in  the  fused  image 
10(c),  Where  the  metric  ignores  the  contribution  of  contrast  in  the  T2  weighted  band  between  gray 
and  white  matter.  As  a  result,  the  gray  to  white  matter  contrast  in  ip(c)  is  solely  dependent  on 
that  in;  9(a,c),  whereas  the  visibility  of  the  tumor  is  due  largely  to  9(b). 


The  contrast  fusion  technique  in  this  paper  can  also  be  used  to  create  color  visualizations 
whose  Utility  surpasses  that  of  current  methods.  Perhaps  the  most  widely  used  and  acknowledged 
procedure  for  producing  color  composites  which  reveal  the  structure  of  multi-band  remotely  sensed 
imagery  is  to  combine  the  top  three  principal  component  images.  Two  standard  techniques  are  to 
use  the  three  top  principal  components  as  the  red,  green  and  blue  channels  of  a  color  image,  or 
alternativelly  as  the  intensity,  hue  and  saturation  components  (in  HSV  coordinates),  in  that  order 
[13].  Figures  11(b)  11(c)  show  examples  of  these  techniques  for  a  4-band  aerial  image  (red,  green, 
blue  arid  near-infrared),  while  11(a)  shows  the  three  bands  in  the  visible  spectrum.  We  can  use  the 
contrast  fusion  method  introduced  in  this  article  to  create  the  intensity  band  of  a  color  composite 
in  the  HSV  coordinate  system.  Figure  11(d)  shows  the  result  of  this  procedure  using  the  hue  and 
saturation  components  from  11(a),  and  11(e)  shows  the  results  of  using  the  hue  and  saturation 
from  11(b).  In  both  cases  we  see  that  the  quality  of  the  resulting  color  composite  is  improved  by 
the  addition  of  the  contrast  fused  image  as  an  intensity  component.  Most  notably,  in  both  cases 
it  becomes  easier  to  discriminate  trees  from  their  shadows  and  their  surrounding  background.  In 
Figure  11(e),  the  image  aquires  a  more  ‘natural’  look,  and  the  road  markings  become  visible. 

In  terms  of  computational  time,  the  algorithm  is  rather  efficient.  For  example,  for  a  4-band 
image  measuring  500  x  500  pixels,  the  time  to  complete  30  iterations  of  (12)  is  just  2.7  seconds 
on  a  Pentium  III  processor  at  700Mhz.  While  our  theoretical  formulation  requires  (12)  to  reach 
a  steady  state,  it  is  usually  the  case  that  visually  acceptable  results  are  reached  quite  quickly, 
and  further  iteration  adds  imperceptible  changes  to  the  fused  image.  Note  that  the  computational 
scheme  proposed  is  of  quadratic  order  in  the  number  of  bands  and  of  linear  order  in  the  number  of 
image  pixels. 


5  Conclusion 

The  need  for  image  fusion  techniques  is  now  greater  than  ever.  Proliferation  of  imaging  modali- 
.ies  and*  sensor  technologies  continues  to  flood  image  analysts  with  increasing  amounts  of  digital 
nformation  for  processing.  New  and  effective  visualization  methods  that  reduce  the  workload  and 
ncrease  the  efficiency  of  analysts  constitute  therefore  a  crucial  area  for  algorithm  development. 

In  thi."!  paper,  we  presented  an  analytic  definition  of  contrast  applicable  to  a  general  multiband 
mage  and  showed  how  our  definition  agrees  with  the  standard  one  in  the  case  of  a  grayscale 
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image.  The  problem  of  grayscale  visualization  with  optimal  contrast  was  stated  then  in  terms 
of  an  image  contrast  vector  Held,  which  allowed  us  to  obtain  a  clean  mathematical  formulation. 
The  mstantiaUon  of  this  formulation  yields  an  algorithm  capable  of  producing  consistently  high 
qu  1  y  visualizations  of  multiband  images.  The  enhanced  feature  content  of  these  visualizations 
IS  a  powerful  aid  for  the  image  analyst,  who  can  more  effectively  select  regions  of  interest  or 
distinctive  image  features  for  further  processing  by  image  understanding  algorithms.  The  range  of 
applicability  extends  from  automatic  target  recognition  and  aerial  mapping  to  detection  of  brain 
esions  in  medical  iinaging.  A  number  of  examples  and  comparisons  support  our  performance  claim 
and  attest  to  the  utility  of  our  formalism.  ^ 


Two  directions  for  further  research  appear  necessary.  First,  it  would  be  interesting  to  solve 
the  reduction  problem  not  for  an  optimal  grayscale  image,  as  we  have  done,  but  for  an  optimal 
projection  ir ;  F”  ^  [0,M].  It  follows  from  Equation  (8)  that  the  Euler-Lagrange  equation  for  this 
problem  is  A(7r  o  c)  =  divV.  Unfortunately  this  problem  is  rather  ill-posed,  so  more  constraints 
should  be  introduced  in  order  to  solve  it  in  a  meaningful  way.  Of  course  by  solving  this  problem 
we  give  up  the  local  consistency  property  of  our  method.  Regardless,  such  a  solution  is  relevant  in 
rases  where  relative  brightness  is  important  in  a  global  sense  and  there  is  little  or  no  need  for  local 
consistency  on  overlapping  sub-images.  This  may  be  the  case  for  current  exploitation  of  multimodal 
nedical  imaging  of  the  head,  where  tissue  attenuation  is  correlated  with  brightness.  Of  course,  the 
rolor  composites  shown  in  Section  4. .3  can  also  be  used  in  this  situation  to  convey  both  contrast 
md  attenuation. 


Secondly,  it  would  be  desirable  to  create  a  hybrid  method  exhibiting  the  best  qualities  of  the 
nultiresolution  techniques  [22,  3,  21,  9],  and  the  method  developed  above.  A  challenging  series 
)f  images  are  the  Ishihara  pseudo-isochromatic  plates  used  for  testing  color  blindness  [17]  (see 
‘igure  12).  We  see  that  in  these  plates,  most  colored  circles  are  immediately  adjacent  to  white 
lackground  pixels,  and  not  to  other  colored  circles.  Thus  the  short-scale  contrast  in  the  image  is  not 
I  lependent  on  the  respective  colors.  The  eye  perceives  the  larger  scale  contrast,  and  so  should  any 
i  ilgorithm  which  seeks  to  extend  the  performance  of  the  human  eye  beyond  the  visible  spectrum. 

I  )f  course,  by  blurring  the  plates  we  can  create  short-scale  contrast  between  the  colored  areas,  and 
his  is  essentially  what  wavelet-based  multiresolution  algorithms  do.  Thus  we  see  that  it  would 
te  beneficial  to  incorporate  the  best  aspects  of  multiscale  techniques  with  those  of  the  algorithm 
i  fttroduced  in  this  article.  Such  a  programme  is  currently  underway,  and  we  hope  to  report  on  the 
]  esults  in  the  near  future. 
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(a)  (b) 


1'Mgure  4;  (a)  Grayscale  version  of  12-band  image  computed  through  PCA.  (b)  Grayscale  version 
if  the  same  image  computed  through  our  algorithm. 
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-band  image  sequence:  (a)  Through  PCA.  (b) 
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(b) 


Jigure  8:  Compaxison  between  (a)  wavelet  fusion  and  (b)  proposed  method,  for  the  fusion  of  two 
ontrast  windows  of  chest  CT. 
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